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ABSTRACT 

Motivated by a class of orbit problems in astrophysics, this paper considers solutions 
to Hill's equation with forcing strength parameters that vary from cycle to cycle. The 
results are generalized to include period variations from cycle to cycle. The development 
of the solutions to the differential equation is governed by a discrete map. For the general 
case of Hill's equation in the unstable limit, we consider separately the case of purely 
positive matrix elements and those with mixed signs; we then find exact expressions, 
bounds, and estimates for the growth rates. We also find exact expressions, estimates, 
and bounds for the infinite products of several 2x2 matrices with random variables in the 
matrix elements. In the limit of sharply spiked forcing terms (the delta function limit), 
we find analytic solutions for each cycle and for the discrete map that matches solutions 
from cycle to cycle; for this case we find the growth rates and the condition for instability 
in the limit of large forcing strength, as well as the widths of the stable/unstable zones. 



1. INTRODUCTION 

This paper presents new results concerning Hill's equation of the form 

-j^ + [Xk + qkQit)]y = 0, (1) 

where the function Q{t) is periodic, so that Q{t + vr) = Q{t), and normalized, so that Qdt = 
1. The parameter is denoted here as the forcing strength, which we consider to be a random 
variable that takes on a new value every cycle (the index k determines the cycle). The parameter 
Afc, which determines the oscillation frequency in the absence of forcing, also varies from cycle to 
cycle. In principal, the duration of the cycle could also vary; our first result (see Theorem 1) shows 
that this generalized case can be reduced to the problem of equation ([T]) . 

Hill's equations [HI] arise in a wide variety of contexts [MW] and hence the consideration of 
random variations in the parameters {qk,^k) is a natural generalization of previous work. This 
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particular form of Hill's equation was motivated by a class of orbit problems in astrophysics [AB]. 
In many astrophysical systems, orbits take place in extended mass distributions with triaxial forms. 
Examples include dark matter halos that envelop galaxies and galaxy clusters, stellar bulges found 
at the centers of spiral galaxies, elliptical galaxies, and young embedded star clusters. These 
systems thus occur over an enormous range of scales, spanning factors of millions in size and 
factors of trillions in mass. Nonetheless, the basic form of the potential is similar [NF, BE, AB] for 
all of these systems, and the corresponding orbit problem represents a sizable fraction of the orbital 
motion that takes place in our universe. In this context, when a test particle (e.g., a star or a dark 
matter particle) orbits within the triaxial potential, motion that is initially confined to a particular 
orbital plane can be unstable to motion in the perpendicular direction [AB]. The equation that 
describes the development of this instability takes the form of equation ([1]). Further, the motion in 
the original orbital plane often displays chaotic behavior, which becomes more extreme as the axis 
ratios of the potential increase [BT]. In this application, the motion in the original orbit plane - in 
particular, the distance to the center of the coordinate system ~ determines the magnitude of the 
forcing strength qk that appears in Hill's equation. The crossing time, which varies from orbit to 
orbit, determines the value of the oscillation parameter Afc. As a result, the chaotic behavior in the 
original orbital plane leads to random forcing effects in the differential equation that determines 
instability of motion out of the plane (see Appendix A for further discussion). 

Given that Hill's equations arise in a wide variety of physical problems [MW], we expect that 
applications with random forcing terms will be common. Although the literature on stochastic 
differential equations is vast (e.g., see the review of [BL]), specific results regarding Hill's equations 
with random forcing terms are relatively rare. 

In this application. Hill's equation is periodic or nearly periodic (we generalize to the case of 
varying periods for the basic cycles), and the forcing strength varies from cycle to cycle. Since 
the forcing strength is fixed over a given cycle, one can solve the Hill's equation for each cycle using 
previously developed methods [MW], and then match the solutions from cycle to cycle using a 
discrete map. As shown below, the long-time solution can be developed by repeated multiplication 
of 2 X 2 matrices that contain a random component in their matrix elements. 

The subject of random matrices, including the long term behavior of their products, is also the 
subject of a great deal of previous work [BL, DE, BD, FK, FU, LR, ME, VI]. In this application, 
however. Hill's equation determines the form of the random matrices and the repeated multiplication 
of this type of matrix represents a new and specific application. Given that instances where analytic 
results can be obtained are relatively rare, this set of solutions adds new examples to the list of 
known cases. 

This paper is organized as follows. In §2, we present the basic formulation of the problem, 
define relevant quantities, and show that aperiodic generalizations of the problem can be reduced to 
random Hill's equations. The following section (§3) presents the main results of the paper: We find 
specific results regarding the growth rates of instability for random Hill's equations in the limit of 
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large forcing strengths (i.e., in the hmit where the equations are robustly unstable). These results 
are presented for purely positive and for mixed signs in the 2x2 matrix map. We also find limiting 
forms and constraints on the growth rates. Finally, we find bounds and estimates for the errors 
incurred by working in the limit of large forcing strengths. This work is related to the general 
existence results of [FU] but provides much more detailed information in our setting. In the next 
section (§4) we consider the limit where the forcing terms are Dirac delta functions; this case allows 
for analytic solutions to the original differential equation. We note that the growth rates calculated 
here (§3) depend on the distribution of the ratios of the principal solutions to equation ([T]), rather 
than (directly) on the distributions of the parameters {Xk,qk)- Using the analytic solutions for the 
delta function limit (§4), we thus gain insight into the transformation between the distributions 
of the input parameters {Xk,qk) and the parameters that specify the growth rates. Finally, we 
conclude, in §5, with a summary and discussion of our results. 



2. FORMULATION 

Definition: A Random Hill's equation is defined here to be of the form given by equation ([T]) 
where the forcing strength and oscillation parameter Afc vary from cycle to cycle. Specifically, 
the parameters qk and \k are stochastic variables that take on new values every cycle < [t] < vr, 
and the values are sampled from known probability distributions Pq{q) and Px{X). 



2.1. Hill's Equation with Fixed Parameters 

Over a single given cycle, a random Hill's equation is equivalent to an ordinary Hill's equation 
and can be solved using known methods [MW] . 

Definition: The growth factor fc per cycle (the Floquet multiplier) is given by the solution to the 
characteristic equation and can be written as 

/.^^±<4^. (2) 

where the discriminant A = A(g, A) is defined by 

A = yi(vr) + ^(7r), (3) 

and where yi and 1/2 are the principal solutions [MW]. 

It follows from Floquet 's Theorem that |A| > 2 is a sufficient condition for instability [MW, 
AS]. In addition, periodic solutions exist when |A| =2. 



- 4 - 



2.2. Random Variations in Forcing Strength 



We now generalize to the case where the forcing strength and oscillation parameter vary 
from cycle to cycle. In other words, we consider each period from t = 0tot = 7rasa cycle, and 
consider the effects of successive cycles with varying values of ((?fc,Afc). 

During any given cycle, the solution can be written as a linear combination of the two principal 
solutions ui and y2- Consider two successive cycles. The first cycle has parameters {qai^a) and 
solution 

fa{t) = aaVlait) + Pay2a{t) , (4) 

where the solutions yia{t) and y2a{t) correspond to those for an ordinary Hill's equation when 
evaluated using the values (ga,Aa). Similarly, for the second cycle with parameters (gb,Ab) the 
solution has the form 

fb{t) = Obyibit) + (3by2b{t) . (5) 

Next we note that the new coefficients and are related to those of the previous cycle through 
the relations 



and I3b = aa—j^{TT) + Pa—7-{-^) 



(6) 



dt ' ' ' " dt 

The new coefficients can thus be considered as a two dimensional vector, and the transformation 
between the coefficients in one cycle and the next cycle is a 2 x 2 matrix. Here we consider the 
case in which the equation is symmetric with respect to the midpoint t = tt/2. This case arises in 
the original orbit problem that motivated this study — the forcing function is determined by the 
orbit as it passes near the center of the potential and this passage is symmetric (or very nearly so) . 
It also makes sense to consider the symmetric case, which is easier, first. Since the Wronskian of 
the original differential equation is unity, the number of independent matrix coefficients is reduced 
further, from four to two. We thus have the following result: 

Proposition 1: The transformation between the coefficients aa,(3a of one cycle and the coefficients 
o^by Pb of the next may be written in the form 

■h ih^-l)/g- 



ab 




Jb. 





9 



h 



aa 

Pa 



M{qa) 



Pa 



(7) 



where the matrix M (defined in the second equality) depends on the values {qa,^a) oind h = yi(7r) 
and g = yi(vr) for a given cycle. 

Proof: This result can be verified by standard matrix multiplication, which yields equation ([6]) 
above. □ 

After A'' cycles with varying values of {qk,^k), the solution retains the general form given 
above, where the coefficients are determined by the product of matrices according to 

N 



OLN 



"0 

Po 



where M^^) = J] Mfc(gfc, A^) • 



(8) 



k=l 
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This formulation thus transforms the original differential equation (with a random element) 
into a discrete map. The properties of the product matrix M^^^ determine whether the solution is 
unstable and the corresponding growth rate. 

Definition: The asymptotic growth rate 700 is that experienced by the system when each cycle 
amplifies the growing solution by the growth factor appropriate for the given value of the forcing 
strength for that cycle, i.e., 



7oo = lim ^- log 



^ 1 / 

n 2^^'^ + 



=1 



(9) 



where = A(gfc, A^) is defined by equation ([3]), and where this expression is evaluated in the limit 
N 00. In this definition, it is understood that if |Ajt| < 2 for a particular cycle, then the growth 
factor is unity for that cycle, resulting in no net contribution to the product (for that cycle). 

Notice that the factor of tt appears in this definition of the growth rate because the original 
Hill's equation is taken to be vr-periodic [MW, AS]. As we show below, the growth rates of the 
differential equation are determined by the growth rates resulting from matrix multiplication. In 
many cases, however, the growth rates for matrix multiplication are given without the factor of vr 
[BL, FK], so there is a mismatch of convention (by a factor of tt) between growth rates of Hill's 
equations and growth rates of matrix multiplication. 



Notice that this expression for the asymptotic growth rate takes the form 

N 



ioo= iim^^^i{qk,h) ^ ii) , (10) 

' k=l 

where 7(^fe,Afc) is the growth rate for a given cycle. The asymptotic growth rate is thus given 
by the expectation value of the growth rate per cycle for a given probability distribution for the 
parameters and A^. 

We note that a given system does not necessarily experience growth at the rate 700 because the 
solutions must remain continuous across the boundaries between subsequent cycles. This require- 
ment implies that the solutions during every cycle will contain an admixture of both the growing 
solution and the decaying solution for that cycle, thereby leading to the possibility of slower growth. 
In some cases, however, the growth rate is larger than 700, i.e., the stochastic component of the 
problem aids and abets the instability. One could also call 'joo the direct growth rate. 



2.3. Generalization to Aperiodic Variations 

Theorem 1: Consider a generalization of Hill's equation so that the cycles are no longer exactly 
TT-periodic. Instead, each cycle has period /ifcvr, where fik is a random variable that averages to 
unity. Then variations in period are equivalent to variations in {q,X), i.e., the problem with three 
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stochastic variables {qk, ^k, l^k) reduces to a t: -periodic problem with only two stochastic variables 



Proof: With this generahzation, the equation of motion takes the form 



dt^ 



+ 



y = o, 



(11) 



where we have normahzed the forcing frequency to have unit amphtude {Q = Q/qk)- Since Q (and 
Q) are vr-periodic, the jth cycle is defined over the time interval < fikt < vr, or < t < ir/fik- We 
can re-scale both the time variable and the "constants" according to 



t fikt , 



^k ^k/ fJ'k — ' 



and 



qk qk/fJ-l = Qj 



so the equation of motion reduces to the familiar form 

d'y 



dt"^ 



+ 



Xj + qjQit) 



y = 0. 



(12) 



(13) 



Thus, the effects of varying period can be incorporated into variations in the forcing strength qk 
and oscillation parameter A^. □ 



3. HILL'S EQUATION IN THE UNSTABLE LIMIT 



In this section we consider Hill's equation in general form (for the delta function limit see §4), 
but restrict our analysis to the case of symmetric potentials so that yi{n) = h = 2/2 (tt)- We also 
consider the highly unstable limit, where we define this limit to correspond to large /i 3> 1. Since 
the 2x2 matrix of the discrete map must have its determinant equal to unity, the matrix of the 
map has the form given by equation ([7]), where the values of h and g depend on the form of the 
forcing potential. 



The discrete map can be rewritten in the general form 



M = h 



1 X 


+ 


"0 


-1/5" 


1/x 1 









In the highly unstable limit /i — > 00, the matrix simplifies to the approximate form 

1/x 1 

where we have defined x = h/g, and where the second equality defines the matrix C. 



(14) 



(15) 



In this problem we are concerned with both the long-time limit N ^ 00, and the "unstable" 
limit h ^ 00. In the first instance considered here we take the unstable limit first, but below we 
analyze precisely the difference between taking the long time limit first and then the unstable limit. 
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3.1. Fixed Matrix of the Discrete Map 

The simplest case occurs when the stochastic component can be separated from the matrix, 
i.e., when the matrix C does not vary from cycle to cycle. This case arises when the Hill's equation 
does not contain a random component; it also arises when the random component can be factored 
out so that X does not vary from cycle to cycle, although the leading factors can vary. In either 
case, the matrix C is fixed. Repeated multiplications of the matrix C are then given by 

C^ = 2^-^C. (16) 

With this result, after cycles the Floquet multiplier (eigenvalue) of the product matrix and the 
corresponding growth rate take the form 

N I ^ 

A = Y\{2hk) and 7 = lim — V log(2/ifc) . (17) 

k=l k=l 

Note that this result applies to the particular case of Hill's equation in the delta function limit (§4), 
where the forcing strength varies from cycle to cycle but the frequency parameter A/; is constant. 
The growth rate in equation (jl7p is equal to the asymptotic growth rate 700 (eq. [9]) for this case. 



3.2. General Results in the Unstable Limit 

We now generalize to the case where the parameters of the differential equation vary from 
cycle to cycle. For a given cycle, the discrete map is specified by a matrix of the form specified by 
equation (fT5l) . where x = = hk/ gk, with varying values from cycle to cycle. The values of Xk 
depend on the parameters {qk,^k) through the original differential equation. After N cycles, the 
product matrix M^'^-' takes the form 



N 



N 



(18) 



k=l 



k=l 



where we have separated out the two parts of the problem. One can show (by induction) that the 
product of matrices Cfc have the form 



N 



c w = n 



k=l 



(19) 



where xi is the value of the variable for the first cycle and where the sums Tix(n) ^_B(7V) ^-re 
given by 

2N-1 2"' ' 



E 



and 



5^ f ■ 

j=l 1 



(20) 
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where the variables rj are ratios of the form 

rj = °^ . (21) 

The ratios rj arise from repeated multiplication of the matrices Ck, and hence the indices lie in 
the range 1 < ai,bi < N . The rj always have the same number of factors in the numerator and the 
denominator, but the number of factors (n) varies from (where rj = 1) up to N/2. This upper 
limit arises because each composite ratio rj has 2n values of xj, which must all be different, and 
because the total number of possible values is A^. 

Next we define a composite variable 

2N-1 

S ^ ^[ST(iV) + ^BiN)] = ^ E (^^- + ^) • (22) 

j=l •? 

With this definition, the (growing) eigenvalue A of the product matrix M^^^ takes the form 

N 



A = Slli2hk) (23) 



k=l 

and the corresponding growth rate of the instability has the form 

N 



7 = lim 



fc=l 



(24) 



The first term is the asymptotic growth rate 700 defined by equation ([9]) and is thus an average of 
the growth rates for the individual cycles. All of the additional information regarding the stochastic 
nature of the differential equation is encapsulated in the second term through the variable S. For 
example, if the composite variable S is finite in the limit A^ 00, then the second term would 
vanish. As shown below, however, the stochastic component can provide a significant contribution 
to the growth rate, and can provide either a stabilizing or destabilizing infiuence. In the limit 
A^ — > 00, we can thus write the growth rate in the from 

7 = 7oo + A7 , (25) 

where we have defined the correction term A7, 

A7= lim -^logS. (26) 



Since the asymptotic growth rate 700 is straightforward to evaluate, the remainder of this 
section focuses on evaluating A7, as well as finding corresponding estimates and constraints. This 
correction term A7 is determined by the discrete map C, whose matrix elements are given by the 
ratios x = h/g, where h and g are determined by the solutions to Hill's equation over one cycle. 
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One should keep in mind that the parameters in the original differential equation are (Xk^Qk)- 
The distribution of these parameters determines the distributions of the principal solutions (the 
distributions of and gk), whereas the distribution of the ratios Xk of these latter quantities 
determines the correction A7 to the growth rate. The problem thus separates into two parts: [1] The 
transformation between the distributions of the parameters {Xk,qk) and the resulting distribution 
of the ratios x^ that define the discrete map, and [2] The growth rate of the discrete map for a 
given distribution of Xk- The following analysis focuses on the latter issue (whereas §4 provides an 
example of the former issue). 



3.3. Growth Rates for Positive Matrix Elements 

This subsection addresses the cases where the ratios x^ that define the discrete map C all 
have the same sign. For this case, the analysis is simplified, and a number of useful results can be 
obtained. 

Theorem 2: Consider the general form of Hill's equation in the unstable limit so that h = yi(7r) = 
2/2 (tt) ^ 1. For the case of positive matrix elements, rj > 0, the growth rate is given by equation 
125\) where the correction term A7 is given by 

1 ^ 12 

where Xji and Xj2 represent two different (independent) samples of the xj variable^ 

Proof: Using the same induction argument that led to equation ()19|) , one finds that from one cycle 
to the next the sums ^^(Ar) and ^b{N) vary according to 

x 

and 

^B{N+l) = ^B{N) + ~^T{N) ■ (29) 

In this notation, the variable x (no subscript) represents the value of the x variable at the current 
cycle, whereas xq represents the value at the initial cycle. The growing eigenvalue of the product 
matrix of equation ([T9|) is given by A = ^^(Af) + '^b{N)- As a result, the eigenvalue (growth factor) 
varies from cycle to cycle according to 



^ ^ (3;/xo)Sg(jv) + (xo/x)Sr(jv) 



(30) 



^Specifically, the index j labels the cycle number, and the indices jl and j2 label two successive samples of the 
X variable; since the stochastic parameters of the differential equations are assumed to be independent from cycle to 
cycle, however, the variables Xji and Xj2 can be any independent samples. 
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The overall growth factor is then determined by the product 



N 

n 



1 + 



(31) 



The growth rate of matrix multiplication is determined by setting the above product equal to 
exp[A^7r7]. The growth rate A7 also includes the factor of 2 per cycle that is included in the 
definition of the asymptotic growth rate 700- We thus find that 



1 



N 



1 + 



{Xji/xj2)^B{N) + {Xj2/Xjl)T.'r(N) 



-'B(N) 



T{N) 



log 2 



vr 



(32) 



Note that this expression provides the correction A7 to the growth rate. The full growth rate is 
given by 7 = 700 + A7 (where 700 is specified by eq. [9] and A7 is specified by eq. [27]). In the 
limit of large A^, the ratio of the sums ^^(Af) ^^id Ti^^^-j approaches unity, almost surely, so that 



^T{N) 



1 



as 



00 . 



(33) 



This result follows from the definition of 'St{n) ^B(Af)- The terms in each of these two sums 
is the product of ratios Xa/xb, and, the terms rj in the first sum T,x(n) the inverse of those 
(1/rj) in the second sum E^^jy)- Since the fundamental variables Xk that make up these ratios, and 
the products of these ratios, are drawn from the same distribution, the above condition (I33p must 
hold. As a consequence, the expression for the growth rate given by equation ([32]) approaches that 
of equation ([27]). □ 

Corollary 2.1: Let be the variance of the composite variable log(xji/xj2) (see Theorem 2). The 
correction to the growth rate is positive semi- definite; specifically, A7 > and A7 ^ in the limit 
(To ^ 0. Further, in the limit of small variance, the growth rate approaches the asymptotic form 
A7^aoV(87r). 

Proof: In the limit of small o"o, we can write Xj = \ + Sj, where \5j\ <C 1. In this limit, equation 
(j27p for the growth rate becomes 



1 ^ 

A7 = hm — V log[2 + 6,1 - 6j2 + S% - 6ji6j2 + 0{6^)] 



log 2 



vr 



(34) 



In the limit \6j\ <C 1, we can expand the logarithm, and the above expression simplifies to the form 



A7 = lim 



N 



Evaluation of the above expression shows that 



1 

A7 = — 

' 27r 



0j2) 



+ 0(5^) 



Svr' 



(35) 



(36) 
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As a result, A7 > 0. In the limit o"o — > 0, all of the Xj approach unity and 5j — > 0; therefore, 
A7 ^ as fJo ^ 0. □ 

Although equation (j27|) is exact, the computation of the expectation value can be difficult in 
practice. As a result, it is useful to have simple constraints on the growth rate in terms of the 
variance of the probability distribution for the variables x^- In particular, a simple bound can be 
derived: 

Theorem 3: Consider the general form of Hill's equation in the unstable limit so that h = yi(7r) = 
2)2 (^) ^ 1- Take the variables rj > 0. Then the growth rate is given by equation Ii25\) and the 
correction term A7 obeys the constraint 

A7 < I , (37) 

where Uq is the variance of the distribution of the variable ^ = \og{x ji / x j2) , and where Xj are 
independent samplings of the ratios Xj = hj/gj. 

Proof: First we define the variable = logr^, where rj is given by equation (j2ip above with 
a fixed value of n. In the limit of large n, the variable has zero mean and will be normally 
distributed. If the variables xj are independent, the variance of the composite variable will be 
given by 

cj| = na'^ . (38) 

As shown below, is order to obtain 2^ terms in the sums ^^^(jv) ^-iid S^j-jy), almost all of the 
variables rj will fall in the large n limit; in addition, n — > 00 in the limit — > 00. As a result, 
we can consider the large n limit to be valid for purposes of evaluating the correction term A7. In 
practice, the variables will not be completely independent, so the actual variance will be smaller 
than that given by equation (I38p : nonetheless, this form can be used to find the desired upper limit. 

Given the large n limit and a log-normal distribution of r^, the expectation values {rj) and 
{1/rj) are given by 

(r,)=exp[naoV2] = (l/r,). (39) 

Note that the variable is normally distributed, and we are taking the expectation value of 
rj = exp^j; since the mean of the exponential is not necessarily equal to the exponential of the 
mean, the above expression contains the (perhaps counterintuitive) factor of 2. As expected, larger 
values of n allow for a wider possible distribution and result in larger expectation values. The 
maximum expectation values thus occur for the largest values of n. Since n < N/2, these results, 
in conjunction with equation (I22p imply that S obeys the constraint 

S < exp[Na^/A] . (40) 

The constraint claimed in equation (I37p then follows immediately. 

Comhinatorics: To complete the argument, we must show that most of the variables rj have 
a large number n of factors (in the limit of large A^). The number of terms in the sums 'St{N) 
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and S^(7v) is large, namely 2^~^. Further, the ratios rj must contain 2n different values of the 
variables x^- The number P{n\N) of different ways to choose the 2n variables for N cycles (and 
hence possible values of Xk) is given by the expression 

Notice that this expression differs from the more familiar binomial coefficient because the values 
of rj depend on whether or not the Xk factors are in the numerator or denominator of the ratio rj. 
Next we note that if n <C A^, then the following chain of inequalities holds for large A^: 

A72n 

^N^^)< (;^«2^-'- (42) 

For large and n <^ N, the central expression increases like a power of A^, whereas the right hand 
expression increases exponentially with A'^. As a result, for n <ti N, there are not enough different 
ways to choose the x^ values to make the required number of composite ratios rj. In order to allow 
for enough different rj, the number n of factors must be large (namely, large enough so that n <^ N 
does not hold) for most of the rj. This conclusion thus justifies our use of the large n limit in the 
proof of Theorem 3 (where we used a log-normal form for the composite distribution to evaluate 
the expectation values {rj) and {l/rj)). □ 

Estimate: Theorem 3 provides an upper bound on the contribution of the correction term A7 
to the overall growth rate. This bound depends on the value of n, which determines the magnitude 
of the expectation value (rj). It is useful to have an estimate of the "typical" size of n. In rough 
terms, the value of n must be large enough so that the number of possible combinations is large 
enough to account for the 2^~^ terms in the sums ST(Af) aiid T,b(n)- For each n, we have P{n\N) 
combinations. As a rough approximation, nP{n\N) accounts for all of the combinations of size less 
than n. If we set nP(n|A^) = 2^, we can solve for the ratio n/N required to have enough terms, 
and find n/N ~ 0.11354 . . . ~ 1/9. As a result, we expect the ratio n/N to lie in the range 

1 n 1 , ^ 

If we use this range of n/N to evaluate the expectation value using equation ([39]) . and estimate 
the growth rate, the upper end of this range provides a rigorous upper bound (Theorem 3). The 
lower end of the range only represents a rough guideline, however, since the variables are not fully 
independent. Nonetheless, it can be used to estimate the expectation values {rj). 

Notice that the upper bound is conservative. Figured] shows a comparison of the actual growth 
rate (from Theorem 2) and the bound (Theorem 3). At large variance, the actual growth rate is 
much less than our bound. In fact, as shown in the following section, in the limit of large variance, 
the growth rate A7 oc do (rather than oc cJq). 

For this numerical experiment, we used a particular form for the x^ variables, namely x^ = 
0.01 + (lOa^fc)", where is a random variable in the range < < 1 and a is a parameter 
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that is chosen to attain varying values of ctq. The exact form of the curve A7((Tq) depends on the 
distribution of the Xk- However, all of the distributions studied result in the general form shown 
in Figure [H and all of the cases show the same agreement between numerical experiments and the 
predictions of Theorem 2. 



3.4. Error Bounds and Estimates 

The analysis presented thus far is valid in the highly unstable limit, as defined at the beginning 
of this section. In other words, we have found an exact solution to the reduced problem, as 
encapsulated in equation (jlSp . In this problem we are taking two limits, the long-time limit N ^ oo 
and the "unstable" limit h ^ oo. In the reduced problem, as analyzed above, we take the limit 
/i — > oo first, and then consider the long-time limit N — > oo. In this subsection, we consider 
the accuracy of this approach by finding bounds (and estimates) for the errors in the growth rates 
incurred from working in the highly unstable limit. In other words, we find bounds on the difference 
between the results for the full problem (with large but finite /i^) and the reduced problem. 

To assess the error budget, we write the general matrix (for the full problem) in the form 

1 X(j) 

1/x 1 

This form is the same as the matrix of the reduced problem (in the unstable limit) except for the 
correction factor (p in the (1,2) matrix element, where </> = (1 — l//i 



M = hB where B 



(44) 
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Let {Aj)b denote the growth rate for the matrix B for the full problem defined in equation 
. Similarly, let {A'y)c denote the growth rate found previously for the reduced problem using 
the matrix C defined in equation (jlSp . Through repeated matrix multiplications, the product of 
matrices B^ will be almost the same as for the product of matrices C^, where the difference is due 
to the continued accumulation of factors (pk- Note that the index k, as introduced here, denotes 
the cycle number, and that all of these quantities vary from cycle to cycle. 

Proposition 2: The error ebc = {^l)c ~ {^i)b introduced by using the reduced form of the 
problem (the matrices Ck) instead of the full problem (the matrices B^^ is bounded by 

<eBC < --(^og(j)k) . (45) 

Proof: Since i;^^ < 1, by definition, we see immediately that the growth rate for the full problem 
is bounded from above by that of the reduced problem, i.e., 

{Aj)b < (A7)c . (46) 
Next we construct a new matrix of the form 

1 x' 

1/x 1 



(pC . (47) 



500 



1000 



1500 



Fig. 1. — Comparison of the bound of Theorem 3 and the prediction of Theorem 2 with results 
from numerical experiments. All cases use matrices of the form given by equation (jlSp . where 
the variables are chosen according to distributions with variance ctq. For each distribution, the 
growth rate A7 due to matrix multiplication is plotted versus the variance of the distribution of the 
composite variable ^ = \og{xj/xk)-, where Xk = VikiT^) / Vikij^) and, similarly, Xj = yij (vr) /yij (vr) . 
The solid curve shows the results obtained by averaging together 1000 realizations of the numerical 
experiments; the overlying dashed curve shows the prediction of Theorem 2. The straight solid line 
shows the upper bound of Theorem 3, i.e., A7 < crQ/(47r). 
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The products of the matrices will be almost the same as those for the matrices B^,, where the 
difference is again due to the inclusion of additional factors of (pk- Since the (pk < 1, we find that 
the growth rate for this benchmark problem is less than (or equal to) that of the full problem, i.e., 
(A7)^ < {A^)b- Further, the growth rate {A'y)A for this new matrix can be found explicitly and 
is given by 



{A'j)a = (A7)c + lim — log 
Af— >oo vriV 



■ N 
.k=l 



1 ^ 

(A7)c + lim — V lo; 
N-^oo vriV ^-^ 



(48) 



k=l 



Combining equations (j46p and (j48p shows that the growth rate for the full problem {A'j)b is 
bounded on both sides and obeys the constraint 



1 



iAj)c + -(log</.fc) < (A7)ij < (A7)c7 . 



vr 



(49) 



Notice that the expectation value {log(j)k) < since (pk < 1- The error ebc introduced by using 
the reduced form of the problem (the matrices C^) instead of the full problem (the matrices B^.) 
is thus bounded by 

(50) 



< EEC < --(^OgCp^ 

vr 



This bound can be made tighter by a factor of 2. Note that the product of two matrices of the 
full problem has the form 



B9B1 



1 + (x2/a;i)(/)2 XiCpi + X2(l)2 
1/Xi + 1/X2 1 + (Xi/X2)</'1 



(51) 



Thus, the product of two matrices contains only linear factors of (pk- As a result, we can define a 
new reference matrix A = that accumulates factors of (pk only half as quickly as the original 

matrix A in the above argument, so that 



AoAi 



1 + X2/X1 Xi + X2 

l/xi + l/3;2 l + Xl/3;2 



^r<prc2c. 



(52) 



The new reference matrix still grows more slowly than the matrix B of the full problem, but the 



product of such matrices accumulates only extra factors of v^^, 
in the above argument results in the tighter bound 



1/2 



Using this reference matrix 



< esc < 



2n 



(log (pk 



(53) 



In the limit where all of the h/. ^ 1, logcpk 
approximate form 



T//i|, and the above bound approaches the 



(54) 



This expression shows that the errors are well controlled. For large but finite h^, the departure 
of the growth rates from those obtained in the highly unstable limit (Theorem 2) are 0(/i^^). □ 
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Given the above considerations, we can write the growth rate (Aj)b for the fuh problem in the 
form 

(A7)b = (A7)c-^(/^,-'), (55) 

where (A7)c is the growth rate for the reduced problem and where K,; is a constant of order unity. 
In the limit of large hk (specifically, for logcpk ~ 1/^1)) the constant is bounded and lies in the 
range < < 1/2. Our numerical exploration of parameter space suggest that Ks ^ 1/4 provides 
a good estimate for the correction term. In any case, however, the correction term depends on hk 
through the quantity (^^^) and decreases with the size of this expectation value. 



3.5. Matrix Elements with Varying Signs 

We now consider the case in which the signs of the variables rj can be either positive or 
negative. Suppose that the system has equal probability of attaining positive and negative factors. 
In the limit N oo, one expects the sums Sj-^jv); ^B(Af) ~^ 0) which would seem to imply no 
growth. However, two effects counteract this tendency. First, the other factor that arises in the 
repeated matrix multiplication diverges in the same limit, i.e., 

N 

JJ(2/ife) ^ oo as iV ^ oo . (56) 

k 

Second, the sums Sy(7v) and Y^b{n) can random walk away from zero with increasing number N of 
cycles, where the effective step length is determined by the variance do defined previously. If the 
random walk is fast enough, the system can be unstable even without considering the diverging 
product of equation (j56p . In order to determine the stability (or instability) of the Hill's equation 
in this case, we must thus determine how the sums S^^^) and S^j-tv) behave with increasing A^. 



Theorem 4: Consider the case of Hill's equation in the unstable limit with both positive and 
negative signs for the matrix elements. Let positive signs occur with probability p and negative signs 
occur with probability 1 — p. Then the general form of the growth rate is given by 



lim — — 



N 



[P' + (1 - Pf] Yl + 1^1) + 2p(l 



N 

P) Y log 

k=l 




log 2 



TT 



(57) 



Proof: The same arguments leading to equation (I32p in the proof of Theorem 2 can be used, 
where the signs of the ratios Xjxjxji must be taken into account. If p is the probability of the Xj 
variables being positive, the probability of the ratio of two variables being positive will be given 
by + (1 — p)^, i.e., the probability of getting either two positive signs or two negative signs. 
The probability of the ratio being negative is then 2p(l — p). With this consideration of signs, the 
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intermediate form of equation (|32p is modified to take the form 



A7 + 



log 2 



Np 



TT 



Ntt 



1 + 



\Xji/Xj2\^B{N) + \Xj2/Xjl\Y:x{N) 
'^B{N) + '^T{N) 

\Xji/Xj2\^B{N) + \Xj2/Xjl\T:x[N) 



(58) 



T(N) 



where Np is the number of terms where the ratios have positive signs and Nq is the number of 
terms where the ratios have negative signs. In the hmit — > 00, we argue (as before) that the sums 
^B{N) ^-iid Sjij-jv) approach the same value. Notice also that the two sums can be either positive or 
negative, but they will both have the same sign (by construction). As a result, we can divide the 
sums out of the expression as before. In the limit N 00, the fraction Np/N p'^ + {1 —p)'^ and 
the fraction Nq/N 2j»(l — p). After some rearrangement, we obtain the form of equation ()57p . 
□ 

Corollary 4.1: Let P{^) denote the probability distribution of the composite variable ^ = Xk/xj, 
and assume that the integral J (i^(dP/(i^) log |^| exists. Then for Hill's equation in the unstable 
limit, and for the case of the variables x^ having mixed signs, in the limit of small variance the 
correction to the growth rate A7 approaches the following limiting form: 



lini A7 = ML^ [log do + Co - log 2] , 

o-Q— >0 vr 

where Cq is a constant that depends on the probability distribution of the variables xt- 



(59) 



Proof: In the limit of small ctq, the variables x^ can be written in the form = 1 + (5^ where 
^ 1. To leading order, the expression of equation (|57p for the growth rate becomes 



A7 + 



log 2 



TT 



lim — — 



N N 

\p' + {p- if] J2 log(2 + Sji - Sj2) + 2p{l -p)Y^ log\dki - 6k2\ 
j=i k=i 



(60) 

In the limit of small variance (Tq — > 0, the variables 6k — > 0, and the above expression reduces to 
the form 

A7 = MlzZ) [(log \5,, - 5k2\) - log2] . 



vr 



We thus need to evaluate the expectation value given by 

(log 14 -5,1) = I d^iogiei 



dP 



(61) 



(62) 



where we have defined the composite variable = Sk — Sj. Notice that in the limit \6\ <^ 1, the 
variance of ^ is cjg. Next we define a dimensionless variable z = so that the integral becomes 



dz^ log{aoz) 



dP 



dz 



dP 



logcTQ / dz— h / dz— logz = logcTo + / dz— log z. (63) 



dz 



dP 



dz 



2 4 6 8 10 

2 
^ 



Fig. 2. — Correction A7 to the growth rate for the case in which the signs of the random variables 
Xk axe both positive and negative. The three curves show the results for a 50/50 distribution 
(bottom), 75/25 (center), and the case of all positive signs (top). For all three cases, the solid 
curves show the results of numerical matrix multiplication, where 1000 realizations of each product 
are averaged together. The overlying dashed curves, which are virtually indistinguishable, show 
the exact results from Theorem 4. 
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Fig. 3. — Convergence of growth rates in the Umit of large variance. The increasing sohd curve 
shows the growth rate as a function of variance for tlie case of all positive signs. The dashed curve 
shows the growth rate for the cased of mixed signs with a 50/50 sign distribution, i.e., p = 1/2. 
The decreasing curve marked by triangles shows the difference between the two curves (where the 
axis on the right applies). 
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As long as the differential probability distribution dP/dz allows the integral in the final expression 
to converge, then I = loguo + Co, where Co is some fixed number that depends only on the shape 
of the probability distribution. This convergence requirement is given by the statement of the 
corollary, so that Corollary 4.1 holds. Notice also that in the limit of small do, the loguo term 
dominates for any fixed Cq, so that A7 ~ 2p(l — p)(log(To)/7r. □ 

Figure [2] shows the growth rates as a function of the variance ctq for the case of mixed signs. 
For the case of positive signs only, p = 1, the correction A7 to the growth rate goes to zero as 
(To ^ 0. For the case of mixed signs, the correction to the growth rate has the form A7 oc logcio as 
implied by Corollary 4.1. 

Sometimes it is useful to explicitly denote when the growth rates under consideration are the 
result of purely positive signs or mixed signs for the variables Xk- Here, we use the notation 
to specify the growth rate when all the signs are positive. Similarly, A71J denotes growth rates for 
the case of mixed signs. 

Corollary 4.2: In the limit of large variance, ctq 00, the growth rates for the case of positive 
signs only and for the case of mixed signs converge, i.e., 

lim A7, = A7p, (64) 

where Ajp denotes the case of all positive signs and A^q denotes the case of mixed signs. 
Proof: The difference in the growth rates for two cases is given by 



2p(l — p) 1 ^ ^ r 1 

A7p-A75 = lim — V log(l + |xji/xj2|) - log 1 - |xji/xj2| , (65) 

vr N^QO I\ — ' L J 



where p is the probability for the sign of x\^ being positive. In the limit of large variance Cq ^ 00, 
the ratios are almost always far from unity. Only the cases with ^ 1 have a 

significant contribution to the sums. For those cases, however, both of the logarithms in the sums 
reduce to the same form, loglxj/x^l, and hence equation ([65|) becomes 

lim A7p - A7,; = — — [(log \xjlxk\) - (log . (66) 



As a result, equation (j64l) is valid. □ 

Corollary 4.3: In the limit of large variance do 00, the growth rate A7 approaches the form 
given by 

lim A7=^Coo, (67) 

(TO— >oo TT 

where Cqo is a constant that depends on the form of the probability distribution for the variables 
Xfc. In general. Coo < 1/2. 

Proof: Let the composite variable ^ = log{xk/xj) have a probability distribution dP/dS,. Since the 
growth rate for the case of mixed signs converges to that for all positive signs in the limit of interest 
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(from Corollary 4.2), we only need to consider the latter case (from Theorem 2). The growth rate 
is then given by the expectation value 

A7 = -/ d^-j-log(l + e^). (68) 

The integral can be separated into the domains ^ < and ^ > 0. For the positive integral, we 
expand the integrand into two terms; for the negative domain, we change variables of integration 
so that ^ — > — ^. We thus obtain the three terms 



A7 = - / di—i + - / di— log 1 + e-«) + - / di— log(l + e-«^ 



(69) 



In the third integral, the probability distribution {dP/d^){(,) = (dP/d^)(— ^); the second and third 
terms will thus be the same since the distribution is symmetric (by construction, the composite 
variable ^ is the difference between two variables logx^ drawn from the same distribution). The 
sum of the second two integrals is bounded from above by log 2 and can be neglected in the limit 
of interest. In the first integral, we change variables according to z = ^/a, so that 

cq dP 
A7^— (z)(^>o) where (^;)(^>o) = / dz—z. (70) 

TT Jo 

Since (1) = 1 and (z^) = 1, by definition, we expect the quantity (^)(5>o) = to be of order 
unity. Further, one can show that Coo as defined here is bounded from above by 1/2. As a result, 
in this limit, we obtain a bound of the form 7r(A7) < ao/2 + log 2. We note that the constant 
Coo cannot be bounded from below (in the absence of further constraints placed on the probability 
distribution dP/d^). □ 

Corollary 4.4: In the limit of large variance Gq S> 1, the difference A(A7) between the growth 
rate for strictly positive signs and that for mixed signs takes the form 

hm A(A7) = Mi^CA, (71) 

(T0-*oo TTCTo 

where Ca is a constant that depends on the form of probability distribution, and where p is the 
probability of positive matrix elements for the case of mixed signs. 

Proof: Using the results from Theorem 2 and Theorem 4 to specify the growth rates for the cases 
of positive signs and mixed signs, respectively, the difference can be written in the form 



A(A7) = ^^^-^ r d^ ^ [log(l + e«) - log |1 - e« 

TT J-00 dt, L 



(72) 



Next we separate the integrals into positive and negative domains and change the integration 
variable for the negative domain (^ — ^). The integral (/) then becomes 



Jo 



, ,l + e-^\ dP fl + e-^\ 
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where P{£,) = P{—^). Since we are working in the large ctq hmit, the variable ^ will be large 
over most of the domain where the integrals have support, so we can expand using e~^ as a small 
parameter. In this case, the integral / becomes 

d^^e-l«l = 2 r dz^e-'^^^^K (74) 

where we have made the substitution z = ^ja. For large <to; the decaying exponential dominates 
the behavior of the integrand. In the limit a — > cxd, the exponential term decays to zero before 
the probability dP/dz changes so that dP/dz — > Ca, where Ca is a constant. The integral thus 
becomes / = 4Ca/o"o, and the difference between the growth rates becomes 

A(A7) = Mi^CA, (75) 

as claimed by Corollary 4.4. □ 

Figure [3] illustrates the behavior implied by the last three Corollaries. In the limit of large 
variance, the growth rates for mixed signs and positive signs only converge (Corollary 4.2). Further, 
growth rates for both cases approach the form A7 oc ctq (as in Corollary 4.3). Finally, the difference 
between the growth rates for the two cases has the characteristic form A(A7) oc l/uo (from 
Corollary 4.4). 

Corollary 4.5: For the case of mixed signs, the crossover point between growing solutions and 
decaying solutions is given by the condition 

[p" + (1 -p)2](log |1 + \xj/xk\\) + 2p(l -p)(log |1 -\xj/xk\\) = log 2. (76) 
Proof: This result follows form Theorem 4 by inspection. □ 

Estimate for the Crossover Condition: Equation (|76|) is difficult to evaluate in practice. In order to 
obtain a rough estimate of the threshold for instability, we can consider the the rj to be independent 
variables and use elementary methods to estimate the conditions necessary for systems with mixed 
signs to be unstable. We first note that the sums ^^(Ar) and ^^(Ar) add up the composite variables 
rj, which are made up of the variables xj (which in turn are set by the form of the original 
differential equation). If the signs of the variables xj are symmetrically distributed, then the signs 
of the composite variables rj are also symmetrically distributed. We can thus focus on the variables 

Since the signs can either be positive or negative, the probability of a net excess of positive (or 
negative) terms is governed by the binomial distribution (which has a gaussian form in the limit of 
large A^) . The probability P of having a net excess of m signs is given by the distribution 

P{m) = exp[-mV2iV5] , (77) 

where Ns is the number of steps in the random walk. The sums Sj>(-Ar) and S^^-^f) have Ns = 2^ 
steps, where is the number of cycles of the Hill's equation. 
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If the net excess of signs of one type is m, the sums are reduced (from those obtained with 
purely positive variables) so that 

S = So^^, (78) 

where Sq is the value of the composite sum obtained when the variables Xj have only one sign. 
The probability of a growing solution is given by 

roo 

Pg= P{m)dm , (79) 

where is the minimum number of steps needed for instability. We can write in the form 

= Nse-^''^^° = exp[iV(log(2) - 7rA7o)] , (80) 

where A70 is the correction to the growth rate for the case of positive signs only. 
The integral can be written in terms of the variable ^ = m/(2A^s)^/^ so that 

2 



Pg 



2 2 

e-'dz, (81) 



where ^ 

= exp[iV(- log 2 - 7rA7o] . (82) 

Thus, the crossover for growth occurs under the condition 

A7o«ilog2/(27r) . (83) 

Keep in mind that this result was derived under the assumption that the variables in the ran- 
dom walk are completely independent. We can derive the above approximate result from a sim- 
pler argument: The sums Sj^^jv) and ^b{n) random walk away from zero according to i\/Ns = 
(r|)^/22A^/2 = exp[ncrg + (iV/2)log2]. As a result, S ^ ex.p[na^ - (iV/2)log2] and hence A7 ^ 
(n/iV)(a2/7r)-(log2)/27r. 



3.6. Specific Results for a Normal Distribution 



In this section we consider the particular case where the composite variable ^ = \og{xk/xj) 
has a normal distribution. Specifically, we let the differential probability distribution take the form 

dP 1_ 



27r(7o 



(84) 



so that (Tg is the variance of the distribution. In order to determine the growth rates, we must 
evaluate the integrals 

1 

(85) 



1 /"^ 

V ^nao J -00 
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In the limit ao 0, the correction part of the growth rate (A7) can be evaluated and has the 

form 



lim A7 = i I [p2 + (1 - pf] ^ + 2p(l - p) [log ao - ^] - 3p{l - p) log 2) , 

(TO— »0 TT O 2 J 



.. _ . - - . . - . (86) 

(70- 

where 7em = 0.577215665 ... is the Euler-Mascheroni constant. Note that for the case of positive 
signs only [p = 1), this expression reduces to the form A7 = c7o/(87r) as in Corollary 2.1. For the 
case of mixed signs, this expression reduces to the form A7 oc logao from Corollary 4.1. 

We can also evaluate the growth rate in the limit of large ctq, and find the asymptotic form 

limA7 = ^^. (87) 
cTo— >oo v27r'^' 



As a result, the constant Coo from Corollary 4.3 is given by Coo = l/\/27r. Note that in this limit, 
the growth rate is independent of the probabilities p and (1 — p) for the variables to have positive 
and negative signs, consistent with Corollary 4.2. In this limit, we can also evaluate the difference 
between the case of positive signs and mixed signs, i.e.. 



Thus, the constant Ca from Corollary 4.4 is given by Ca = l/v27r for the case of a normal 
distribution. Note that although Ca = Coo for this particular example, these constants will not be 
the same in general. 

Finally, for the case of purely positive signs, we can connect the limiting forms for small 
variance and large variance to construct a rough approximation for the whole range of ao, i.e., 

A7 - . (89) 

8 + ^/2^0-0 

This simple expression, which is exact in the limits (Tq — > and ao 00, has a maximum error of 
about 18% over the entire range of ao- 



3.7. Matrix Decomposition for Small Variance 

For completeness, and as a consistency check, we can study the growth rates by breaking the 
transformation matrix into separate parts. In this section we consider the case of small variance 
(see Appendix B for an alternate, more general, separation). In the limit of small variance, ctq <C 1, 
the variables Xk only have small departures from unity and can be written in the form 

Xk = l + 5k, (90) 

where \5k\ <C 1. The matrices of the discrete map can then be decomposed into two parts so that 



(91) 
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where = ±1 is the sign of the A;th term, and where 

'1 Sk 
_Sk 1 . 

The matrices A^. and B^. have simple multiphcative properties. In particular, 



and 



1 
-1 



2Ao if Sj 



and 



Sk , 



B? 



but 



-Bi. 



and 



if Sj / Sk , 



I. 



(92) 



(93) 



(94) 



'k -"-"A; ) aiiu jj^ 

The product matrix Y\ Ck will contain long strings of matrices A^ and B^ multiplied by each other. 
If any two matrices A^, have opposite signs in such a multiplication string, then the product of the 
two matrices will be zero and the entire string will vanish. As a result, after a large number N 
of cycles, the only matrices that are guaranteed to survive in the product are those with only B^ 
factors and those with only one A^ factor. Although it is possible for strings with larger numbers of 
Afc to survive, it becomes increasingly unlikely (exponentially) as the number of factors increases. 
To a good approximation, the eigenvalue of the resulting product matrix will be given by the 
product 



N 



A(^) - n • 



(95) 



k=l 



We could correct for the possibility of longer surviving strings of A^ by multiplying by a factor of 
order unity; however, such a factor would have a vanishing contribution to the growth rate. The 
corresponding growth rate thus takes the form 

2p{l-p) 



A7 



lim 

N^oo 



Ntt 



(96) 



k=l 



where the factor 2p{l — p) arises because the matrices with all positive signs lead to a zero growth 
rate in the limit (Tq 0, so only the fraction of the cases with mixed signs contribute. Next we 
note that the sum converges to an expectation value 

f dP 

{M) = J dS—log\6\. (97) 

Next we make the substitution z = 6/ao and rewrite the integral in the form 

f dP f dP 
{M) = (yo J dz—+ / dz— log z. (98) 

0, the first term dominates and the growth rate (to leading order) 

2p{l-p) 



In the limit of interest, do 
approaches the form 



A7 



vr 



log do . 



(99) 



This form agrees with the leading order expression found earlier in Corollary 4.1 (see also Figure 
[21 which shows the growth rate as a function of the variance) . 
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4. HILL'S EQUATION IN THE DELTA FUNCTION LIMIT 

In many physical applications, including the astrophysical orbit problem that motivated this 
analysis, we can consider the forcing potential to be sufficiently sharp so that Q(t) can be considered 
as a Dirac delta function. For this limit, we specify the main equation considered in this section: 

Definition: Hill's equation in the delta function limit is defined to have the form 

^ + [X + q5{[t]-7T/2)]y = 0, (100) 

where q measures the strength of the forcing potential and where 6{t) is the Dirac delta function. 
In this form, the time variable is scaled so that the period of one cycle is vr. The argument of the 
delta function is written in terms of [t], which corresponds to the time variable mod-7r, so that the 
forcing potential is 7r-periodic. 

This form of Hill's equation allows for analytic solutions, as outlined below, which can be used 
to further elucidate the instability for random Hill's equations. In particular, in this case, we can 
solve for the transformation between the variables {Xk,Qk) that appear in Hill's equation and the 
derived composite variables that determine the growth rates. 



4.1. Principal Solutions 

To start the analysis, we first construct the principal solutions to equation (jlOOp for a particular 
cycle with given values of forcing strength q and oscillation parameter A. The equation has two 
linearly independent solutions yi{t) and y2{t)-, which are defined through their initial conditions 

yi(0) = l, ^(0) = 0, and yaW = 0, ^(0) = 1. (101) 
at at 

The first solution yi has the generic form 

yi{t) = cos \/~\t for < t < 7r/2 , (102) 

and 

yi{t) = A cos y/\t + Bsni\/\t for 7r/2<t<7r, (103) 

where A and B are constants that are determined by matching the solutions across the delta 
function at t = 7r/2. We define 9 = \/A7r/2 and find 

A = 1 + (g/VA)sin6lcos6l and B = -{q/ \/\) cos^ 9 . (104) 

Similarly, the second solution y2 has the form 



y2{t) = sin\/At for < t < 7r/2, (105) 
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and 



where 



2/2 



{t) = CcosVXt + Dsin^/Xt for 7r/2 < f < 



C={q/X)sm'^9 and 1) = ^ - (g/A) sin 6* cos ( 

V A 



(106) 
(107) 



For the case of constant parameters {q,X), we can find the criterion for instabihty and the 
growth rate for unstable solutions. Since the forcing potential is symmetric, yi{ir) = dy2/dt{TT), 
from Theorem 1.1 of [MW]. The resulting criterion for instability reduces to the form 



H 



and the growth rate 7 is given by 



2VA 



sm 



(■v/Att) — cos(\/A- 



> 1, 



7 



vr 



(108) 



(109) 



In the delta function limit, the solution to Hill's equation is thus specified by two parameters: the 
frequency parameter A and the forcing strength q. Figure H] shows the plane of possible parameter 
space for Hill's equation in this limit, with the unstable regions shaded. Note that a large fraction 
of the plane is unstable. 



4.2. Random Variations in the Forcing Strength 

We now generalize to the case where the forcing strength q varies from cycle to cycle, but the 
oscillation parameter A is fixed. This version of the problem describes orbits in triaxial, extended 
mass distributions [AB] and is thus of interest in astrophysics. As outlined in §2.2, the solutions 
from cycle to cycle are connected by the transformation matrix given by equation ([7|). Here, the 
matrix elements are given by 

/i = cos(\/Avr) ^ sin(\/A'/r) and 5 = — \/Asin(-v/A7r) — gcos^(\/A7r/2) . (HO) 

2v A 

Theorem 5: Consider a random Hill's equation in the delta function limit. For the case of fixed 
X, the growth rate of instability approaches the asymptotic growth rate 700 in the highly unstable 
limit q/^/X ^ 1, where the correction term has the following order: 




Corollary 5.1: In the delta function limit, the random Hill's equation with fixed X is unstable 
when the asymptotic growth rate joo > 0. 
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Fig. 4. — Regions of instability for Hill's equation in the delta function limit. The shaded regions 
show the values of {X,q) that correspond to exponentially growing (unstable) solutions, which 
represent unstable growth of the perpendicular coordinate for orbits in our triaxial potential that 
are initial confined to one of the principal planes. 
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Remark 5.2: Note that j^o > requires only that a non-vanishing fraction of the cycles are 
unstable. 

Proof: For this version of the problem, the matrix M represents the transition from one cycle to 
the next, where the solutions are written as linear combinations of yi and y2 for the given cycle. In 
other words, this transformation operates in the (2/1,2/2) basis of solutions. However, one can also 
consider the purely growing and decaying solutions, which we denote here as /+ and /_. 

For a given cycle, the eigenvectors V± of the matrix M take the form 

V± 



1 

±g/k 



(112) 



where the +(— ) sign refers to the growing (decaying) solution. The eigenvalues have the form A± 
= hzizk, where k = (/i^ — 1)^^"^. Keep in mind that h = yi(7r) and g = y2{'^), and that A_ = l/A^-. 
We can write any general solution in the form 



f = AV+ + BV^ , 



(113) 



where the coefficients [A, B) are related to the coefficients (a, (3) in the first basis through the 
transformation 



(114) 



^1 _ 1 [1 k/g 1 \a 
b\ ~ 2 [1 -k/g\ [15 

In the basis of eigenvectors, the action of the differential equation over any cycle is to amplify 
growing solution (eigenvector) and attenuate the decaying solution, and this action can be written 
as the matrix transformation 

\ An r A , n 1 r 4 1 

(115) 



A' 
B' 



A 
B 



A+ 
A_ 

At the end of the cycle, we can transform back to the original basis through the inverse of the 
transformation (|114p . As a result, the original matrix M can be decomposed into three components 
so that 

iri 1 iTA, nin hin ^ 

(116) 



1 


1 


1 




"A+ ■ 




"1 k/g - 


2 


.g/k 


-9/k_ 




A„ 




I -k/g_ 



For each cycle, the values of (g. A) can vary. The next cycle will have a new matrix of the same 
general form, with the matrix elements specified by ((/', A'). 

We now shift our view to the basis of eigenvectors, so that each cycle amplifies the growing 
solution. Between the applications of the amplification factors, the action of successive cycles 
"rotates" the solution according to a transition matrix of the form 



T(g,A;g',A') 



k'/g' 
-k'/g' 



1 

g/k 



1 

-g/k 



i + n i-iz 
i-n i + n 



[117) 



where the primes denote the second cycle and where we have defined TZ = k'g/{kg'). For the case 
in which successive cycles have the same values of the original parameters (g. A), the transition 
matrix T becomes the identity matrix (as expected). 
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For simplicity, we now specialize to the case where A is held constant from cycle to cycle, but 
the forcing strength q varies. We can evaluate the transition matrix for the case in which Hill's 
equation hes in the delta function limit and where we also take the limit ql\f\ ^ 1. In this regime. 



7^= 1 + 



q-(i 2VA 1-2 cos( VAtt) 



sin 



(VAtt) 



+ 



A 



= 1 + 2(5. 



(118) 



Note that 7?. = 1 + 2(5 to leading order, where (5 (defined through the above relation) is small 
compared to unity and the sign of 8 can be both positive and negative. Thus, not only is the 
parameter b small, but it can average to zero. Repeated iterations of the mapping lead to the (1,1) 
matrix element growing according to the product 



AT 



(1,1) 



fe=i 



N 



.fe=i 



1+ 



N 
fe=l 



N 



5k + Y,0{5l) 



fe=i 



(119) 



The other matrix elements are of lower order (in powers of l/g) so that to leading order the growing 
eigenvalue of the product matrix is equal to the (1,1) matrix element. Further, for sufHciently well- 
behaved distributions of the parameter g, the sum of 6k averages to zero as iV — > oo. The growth 
rate is thus given by 



N 



N 



) = 7o 



fe=l k=l ^ 

The condition required for the 5k to average to zero can be expressed in the form 



(120) 



1™ 



1 ^ / 1 1 

qq' \' 



0. 



(121) 



k=l 



which will hold provided that the expectation value (l/q) exists. This constraint is nontrivial, 
in that a uniform probability distribution P{q) = constant that extends to g = will produce 
a divergent expectation value for (l/q). Fortunately, in the physical application that motivated 
this analysis, the value of q is determined by the distance to the center of an orbit (appropriately 
weighted) so that the minimum value of q corresponds to the maximum value of the distance. Since 
physical orbits have a maximum outer turning point (due to conservation of energy), physical orbit 
problems will satisfy the required constraint on the probability distribution. □ 



4.3. Second Matrix Decomposition 



Another way to decompose the transformation matrix is to separate it into two separate rota- 
tions, one part that is independent of the forcing strength q, and another that is proportional to q. 
We can thus write the matrix in the form 



M(g,A) 



2VX 



B = 



cos 2^ (sin2^)/VA 
-VAsin26' cos 26' 



sin 261 



2y/X [2VA 008^6' 



sin 26* 



(122) 
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where the second equahty defines the matrices A and B. With these definitions, one finds that 

A^{9) = A{Ne) and (6) = {2sm2e f-^B{9) , (123) 

where we again take A to be constant from cycle to cycle. As a result, after N cycles, the effective 
transformation matrix can be written in the form 



N 

MW = n(A--^B) 

k=l 



(124) 



(2sin20)^-iB(e). 



In the asymptotic limit g/\/A — > oo, the matrix approaches the form 

" N 

mW=(-i)^ n^TTi 

.k=i A_ 

The condition for stability takes the form |TrM^^-*| > 2, i.e. 



(125) 



■ N 

Hqk 

.k=l . 



sin 26* 

~7r 



N 



> 1. 



(126) 



When the system is unstable, the factor on the left hand side of this equation represents the growth 
factor over the entire set of N cycles. The growth rate 7 is thus given by 



7 = lim loe 



N 



n « 



.fc=i 



sin 26* 
TT 



N 



, 1 / sin 26* 

lim > log Ol. r=r- 



k=l 



(127) 



Since = (/^(sin 20)/\/A in this asymptotic limit, the above expression for the growth rate can 
be rewritten in the form 



1 ^ 1 ^ 



k=l 



N^oo N 



(128) 



k=l 



in agreement with Theorem 5. 



4.4. Width of Stable and Unstable Zones 

In the plane of parameters (e.g., Figure (H), the width of the stable and unstable zones can be 
found for the delta function limit. In this case, the leading edge of the zone of stability is given by 
the condition 

e = VXtt = nvr , (129) 

where n is an integer that can be used to label the zone in question. The beginning of the next 
unstable zone is given by the condition \h\ = 1. In the limit of large q ^ 1, the width of the stable 
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regime is narrow, and the boundary will fall at 6 = nir + ip, where ip is small. In particular, ip will 
be smaller than tt/2, so that the angle 9 will lie in either the first or third quadrant, which in turn 
implies that s'mO and cos6 have the same sign. As a result, the condition at the boundary takes 
the form 

g ^l + cos^^2^ (130) 
2VA smip if 

If we solve this expression for ip and use the definition ip = 9 — nir, we can solve for the value of A 
at the boundary of the zone, i.e.. 



n 



2 



A ~ — ; — TTT ~ n 



(1 - 4M)2 

The width of the stable zone can then be expressed in the form 

n2 



1 + — + 0(g-2) . (131) 
qTT 



AA = — . (132) 
irq 

For any finite q, there exists a zone number n such that n? > q and the width of the zone becomes 
wide. In the limit g — > oo, the zones are narrow for all finite n. 

Note that when the forcing strength qj. varies from cycle to cycle, we can define the expectation 
value of the zone widths, 

8n2 / 1 



(AA) = — (-). (133) 
\qk/ 

This expectation value exists under the same conditions required for Theorem 5 to be valid. 



4.5. Variations in (Xk^qk) and Connection to the General Case 

As outlined earlier, the growth rates A7 depend on the ratios of the principal solutions, rather 
than on the input parameters (A^, q^) that appear in the original differential equation ([1]). Since we 
have analytic expressions for the principal solutions in the delta function limit, we can study the 
relationship between the distributions of the fundamental parameters {\k,qk) and the distribution 
of the composite variable ^ = log{xk/xj) that appears in the Theorems of this paper. 

As a starting point, we first consider the limiting case where q^ ^ 00 and the parameter A^ 
is allowed to vary. We also focus the discussion on the correction A7 to the growth rate, which 
depends on the ratios x^- In this limit, using equation (jllOp . we see the variables Xk reduce to the 
simple form 

TT srn^fc 
9k 1 + cos 9k 

where 9k = ^/Xtt. In this case the distribution of ^ = log{xk/xj) depends only on the distribution 
of the angles 9^, which is equivalent to the distribution of A^. Since the Xj and x^ are drawn 
independently from the same distribution (of 9k), the variance of the composite variable ctq = 2c7^, 
where cj^ is the variance of log Xk ■ 
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As a benchmark case, we consider the distribution of 9 to be uniformly distributed over the 
interval [0,27r]. For this example, 



r^^ de 


log 




lo 2^ 





sml 



1 + cos ( 



- 2 


" /"^^ d9 , / 











vr sin( 



1 + cos ( 



(135) 



Numerical evaluation indicates that o"o ~ 2.159. Further, the correction to the growth rate is 
bounded by A7 < Uq/^Att) ~ 0.371 and is expected to be given approximately by A7 ~ 0.13. In 
this limit we expect the asymptotic growth rate to dominate. For example, if qk ~ 1000, a typical 
value for one class of astrophysical orbits [AB], then 700 ~ 2, which is an order of magnitude greater 
than A7. Note that in the limit of large (but finite) g^, the corrections to equation (jl34p are of 
order 0{l/qk), which will be small, so that the variance ctq of the composite variable ^ will be 
nearly independent of the distribution of Qk in this limit. 

As another way to illustrate the transformation between the (A^, qk) and the matrix elements 
Xk, we consider the case of fixed and large but finite (and varying) values of qk- We are thus 
confining the parameter space in Figure H] to a particular vertical line, which is chosen to be in an 
unstable band. We thus define 6 = ^/Xtt, and the Xk take the form 



Xk 



gfc(7r/0) sinO — 2 cos 6 
qkil + cos 9) /2+ {6 /it) sine ' 



(136) 



For purposes of illustration, we can make a further simplification by taking 9 to have a particular 
value; for example, if = 7r/2, the Xk are given by 



Xk 



For this case, the relevant composite variable ^ is given by 

qj + l 



i = log 



qj Qk + 1_ 



(137) 



(138) 



where qj and qk are the values for two successive cycles. In the limit of large qj,qk 3> 1, the 
composite variable takes the approximate form ^ ~ {Qk~Qj) I {QkQj) which illustrates the relationship 
between the original variables (only the qk in this example) and the Xk, or the composite variable 
^, that appear in the growth rates. 

Before leaving this section, we note that the more general case of Hill's equation with a square 
barrier of finite width can also be solved analytically (e.g., let Q{t) = 1/w for a finite time interval 
of width At = w, with Q{t) = otherwise). For this case, in the limit of large qk, the solution for 
hk takes the form 

^l/2 ( Ik 



\hk\ oc sm{wqk 



wXi 



1/2 



(139) 



In the limit of large but finite qk and vanishing width u; — > 0, we recover the result from the delta 
function limit, i.e., the dependence on the width w drops out and \hk\ oc qk- In the limit of finite w 
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and large [specifically, when {wq^) <^ 1 does not hold], then oc ^/qk- This example vindicates 
our expectation that large q^ should lead to large /i^, but the dependence depends on the shape of 
the barrier Q{t). An interesting problem for further study is to place constraints on the behavior 
of the matrix elements hk (and g^) as a function of the forcing strengths q^ for general Q{t). 

5. DISCUSSION AND CONCLUSION 

This paper has considered Hill's equation with forcing strengths and oscillation parameters 
that vary from cycle to cycle. We denote such cases as random Hill's equations. Our first result 
is that Hill-like equations where the period is not constant, but rather varies from cycle to cycle, 
can be reduced to a random Hill's equation (Theorem 1). The rest of the paper thus focuses on 
random Hill's equations, specifically, general equations in the unstable limit (§3) and the particular 
cases of the delta function limit (§4), where the solutions can be determined in terms of elementary 
functions. 

For a general Hill's equation in the limit of a large forcing parameter, we have found general 
results governing instability. In all cases, the growth rates depend on the distribution of values 
for the elements of the transition matrix that maps the solution for one cycle onto the next. 
The relevant composite variable ^ is determined by the principal solutions via the relation ^ = 
log[yiA;(7r)yij(7r)/yifc(7r)yij(7r)], where k and j denote two successive cycles; our results are then 
presented in terms of the variance ctq of the distribution of ^. The growth rate can be separated into 
two parts, the asymptotic growth rate 700 that would result if each cycle grew at the rate appropriate 
for an ordinary Hill's equation, and the correction term A7 that results from matching the solutions 
across cycles. The asymptotic growth rate 700 is determined by the appropriate average of the 
growth rates for individual cycles (see eqs. [9] and |10j). In contrast, the correction term A7 
results from a type of random walk behavior and depends on the variance of the distribution of the 
composite variable ^ defined above. 

For the case of purely positive matrix elements, the correction term A7 has a simple form 
(Theorem 2), and is positive semi-definite and bounded from above and below. In the limit of 
small variance, the correction term A7 oc ctq, whereas in the limit of large variance, A7 oc ctq. For 
all ctQ) the correction term to the growth rate is bounded by A7 < (Jq/4:1t (Theorem 3). A sharper 
bound could be obtained in the future. 

For the case of matrix elements with varying signs, we have found the growth rate of instability 
(Theorem 4) , where the results depend on the probability p of the matrix elements having a positive 
sign. In the limit of small variance, the correction term A7 is always negative and approaches the 
form A7 oc log (To (unless p = 1, where A7 ^ in this limit). As a result, the total growth rate 
7 = 7oo + A7 will always be negative - and hence the system will be stable - for sufficiently small 
variance ctq and any admixture of mixed signs. In the opposite limit of large variance, the growth 
rate for mixed signs and that for purely positive signs converge, with both cases approaching the 



- 35 - 



form A7 oc (t; the difference between the growth rates for the two cases decreases as A(A7) oc I/ctq. 

For the delta function hmit, we can find the solution explicitly for each cycle, and thus an- 
alytically define the matrix elements of the discrete map that develops the solution (eq. [7] and 
|110j ). For the case in which only the forcing strength varies, the growth rate of the general solution 
approaches the asymptotic growth rate (eq. [9]), which represents the growth the solution would 
have if every cycle grows at the rate appropriate for a standard (non-stochastic) Hill's equation. We 
have calculated the widths of the stable and unstable zones for Hill's equation in the limit of delta 
function forcing and large growth rates, which represents a specific case of the results presented in 
[WK], where this specific case includes random forcing terms. Finally, we have used the analytic 
solutions for the delta function limit to illustrate the transformation between the original variables 
i^k^Qk) that appear in Hill's equation and the variables Xk that determine the growth rates (§4.5). 

Although this paper takes a step forward in our understanding of Hill's equation (in particular, 
generalizing it to include random forcing terms) and the multiplication of random matrices (of the 
particular form motivated by Hill's equation), additional work along these lines can be carried out. 
The analysis presented herein works primarily in the limit of large q^, where the solutions are highly 
unstable, although we have bounded the errors incurred by working in this limit. Nonetheless, the 
case in which some cycles have stable solutions, while others have unstable solutions, should be 
considered in greater detail. This paper presents bounds on the correction term A7 to the growth 
rate, but a sharper bound could be found. In the treatment of this paper, we considered the 
probability distribution of the composite variable ^ = log{xk/xj) to be symmetric, which implies 
that Xk and xj are independently drawn from their distribution. In future work, correlations 
between successive cycles can be considered and would lead to asymmetric probability distributions. 
Most of the results of this paper are presented in terms of the distributions of the composite variables 
Xk, rather than the original parameters that appear in Hill's equation; the transformation between 
the distributions of the {Xk,qk) and the Xk thus represents another interesting problem for future 
study. Another case of interest we intend to consider is the case where Q{t) takes the form a 
finite Fourier series. Finally, the relationship between solutions to random Hill's equations and the 
multiplication of random matrices should be explored in greater generality. 

Random Hill's equations, and the properties of their solutions, have a wide variety of appli- 
cations. The original motivation for this work was a class of orbit problems in astrophysics. In 
that context, many astrophysical systems — young embedded star clusters, galactic bulges, and 
dark matter halos — are essentially triaxial extended mass distributions. Orbits within these mass 
distributions are often chaotic; further, when motion is initially confined to a plane, the equation 
of motion for the perpendicular direction is described by a random Hill's equation. The instability 
explored here thus determines how quickly an orbiting body will explore the perpendicular direc- 
tion. For example, this class of behavior occurs in young embedded star clusters, which begin in 
highly flattened configurations but quickly become rounder, in part due to the instability described 
here. Dark matter halos are found (numerically) to display nearly universal forms for their density 
distributions [NF, BE], but an a priori explanation for this form remains lacking. Since the orbits 
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of dark matter particles will be subject to the instability studied herein, random Hill's equations 

must play a role in the explanation. As yet another example, galactic bulges often harbor super- 
massive black holes at their centers; the resulting stellar orbits, including the instability considered 
here, play a role in feeding stars into the central black hole. Finally, we note that in addition to 
astrophysical applications, random Hill's equations are likely to arise in a number of other settings. 
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Appendix A: Astrophysical Motivation 



This Appendix outlines the original astrophysical problem that motivated this study of Hill's 
equation with random forcing. In the the initial setting, the goal was to understand orbits in 
potentials resulting from a density profile of the form 



f{m) 



m 



(Al) 



where po is a density scale. This form arises in many different astrophysical contexts, including 
dark matter halos, galactic bulges, and young embedded star clusters. The density field is constant 
on ellipsoids and the variable m has a triaxial form 



2 

0^ 



(A2) 



where, without loss of generality, a > 6 > c > 0. The radial coordinate ^ is given by = x^+y^+z^. 
The function f{m) is assumed to approach unity as m ^ so that the density profile approaches 
the form p ~ 1/m. For this inner limit, one can find an analytic form for both the potential and 
the force terms [AB]. For purposes of illustration, we write the force terms for the three spatial 
directions in the form 



T 2y 

' \m\ 



2x 



In 



2F(o)VT+2T - Aa' 



a2[2F(a)C + A-2a2^2] 



sm 



. 2r/b2-A ^ 

sm , = 



(A3) 
(A4) 
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2F(c)A/r + 2r - k(? 



(A5) 



c2[2F(c)e + A-2c2e2] 

The coefficients in the numerators are given by the following quadratic functions of the coordinates: 

A^(62 + ^2)^2^(^2^^2y ^(^2^^2)^2 p ^ ^2^2^2 ^ ^2^2y2 ^ ^2j2^2 ^ (^g) 

and the remaining function F is defined by 

F[a) = {i^a^ - ko? + r] . (A7) 



Equations (|A3I - IA7h define the force terms that determine the orbital motion of a test particle 
moving in the potential under consideration (i.e., that resulting from a triaxial density distribution 
of the form |Alj ). The work of [AB] shows that when the orbit begins in any of the three principal 
planes, the motion is (usually) highly unstable to perturbations in the perpendicular direction. For 
example, for an orbit initially confined to the x — z plane, the amplitude of the y coordinate will 
(usually) grow exponentially with time. In the limit of small y, the equation of motion for the 
perpendicular coordinate simplifies to the form 

(Py 2 ,2 4/6 , , , 

-n2+^yy = ^ '^liere iVy = —====== — . (A8) 

Here, the time evolution of the coordinates (x, z) is determined by the orbit in the original x — z 
plane. Since the orbital motion is nearly periodic, the (x, z) dependence of u)y represents a periodic 
forcing term. The forcing strengths, and hence the parameters qk appearing in Hill's equation ([T]), 
are thus determined by the inner turning points of the orbit (with appropriate weighting from the 
axis parameters [a, b,c\). Further, since the orbit in the initial plane often exhibits chaotic behavior, 
the distance of closest approach of the orbit, and hence the strength of the forcing, varies from 
cycle to cycle. The orbit also has outer turning points, which provide a minimum value of Wy, 
which defines the unforced oscillation frequency appearing in Hill's equation ([1]). As a result, 
the equation of motion (lASP for the perpendicular coordinate has the form of Hill's equation, where 
the period, the forcing strength, and the oscillation frequency generally vary from cycle to cycle. 



Appendix B: Growth Rate for an Ancillary Matrix 

In this Appendix, we separate the transformation matrix for the general case (not in the limit 
of small variance) and find the growth rate for one of the matrices. We include this result because 
examples where one can explicitly find the growth rates (Lyapunov exponents) for random matrices 
are rare. Specifically, the transition matrix can be written in the form given by equation (j9ip . where 
the second term in the sum has the form 



Skixk - l)Bfc where 



1 

-1/xk 



(Bl) 
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Note that any pair of matrices with opposite signs will vanish, and so will all subsequent 
products. 

The products of the second term (the matrices along with the leading factor) have a well-defined 
growth rate: 

Proposition 3: The growth rate of matrix multiplication for the matrix = {xk — l)Bjfc is given 
by 



7r = lim 

N^oo 2-kN 



1 f ^ iV ^ 

— I J^log|xfc- 1| + J];iog|l/xj - 1| \ 

[fc=l j=l J 



(B2) 



Proof: The products of the matrices B^ follow cycles as shown by the first three nontrivial cases: 



BoB 





-1/X2 



B3B2B1 



-1/X2 

l/{xix^) 



and 



B4B3B2B1 



(B3) 
(B4) 



l/(a;iX3) 
1/(2:23:4) 

Thus, the even products of the matrices are diagonal matrices, whereas the odd products produce 
matrices with only off-diagonal elements. As a result, the product matrix will approach the form 



N 



\fc=i 

where we have defined 



-fodd 

n p 

^ ^ even 



N 



or 



odd 



N 



\k=l 



N 





-fodd 



-Pr. 



even 





and Pp, 



k=l.odd 



n i 



(B5) 
(B6) 



k=2,even 



For A'" even (odd), the eigenvalues are A = -Peven^-Podd (A = ±z-v/-feven-Podd)- Since |-Peven| = |-Podd| 
in the limit — > 00, the eigenvalues (and hence the growth rates) have the same magnitudes in 
either case. To compute the growth rate jb, we need to account for the fact that only half of 
the factors (either the even or odd terms) appear in the products -Podd and Peven- After some 
rearrangement, we obtain equation ()B2p . □ 
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